This analysis examines the effects of wounding on respiration rates in two coral species: - Porites sp. (por) - Acropora pulchra (acr)
Experimental Design: - 36 total corals (18 per species) - 3 treatments per species (n=6 per treatment): - Control (0): No wounding - Wound Type 1 (1): First wound method - Wound Type 2 (2): Second wound method - Wound date: 2023-05-27 - Measurements over 4 timepoints spanning 23 days
Physiological Measurements: 1. Respiration rates (O2 consumption) 2. Photosynthesis rates (O2 production) 3. Growth (buoyant weight) 4. Photosynthetic efficiency (Fv/Fm via PAM fluorometry) 5. Surface area (wax dipping method)
suppressPackageStartupMessages({
library(tidyverse)
library(lmerTest)
library(emmeans)
library(rstatix)
library(janitor)
library(ggthemes)
library(knitr)
library(kableExtra)
library(pbkrtest)
library(sjPlot)
library(cowplot)
library(patchwork)
library(here)
})
# Source custom functions (loads remaining packages without conflicts)
source("scripts/analysis_functions.R")
# Set theme for all plots
theme_set(theme_few(base_size = 12))
# Define consistent color scheme and labels for all plots
treatment_colors <- c("0" = "#2E86AB", "1" = "#A23B72", "2" = "#F18F01")
treatment_labels <- c("0" = "Control", "1" = "Small Wound", "2" = "Large Wound")
# Helper functions that may need Rmisc
if (!requireNamespace("Rmisc", quietly = TRUE)) {
install.packages("Rmisc")
}# Load sample metadata
sample_info <- read_csv("data/metadata/sample_info.csv", show_col_types = FALSE)
# Convert to factors
sample_info <- sample_info %>%
mutate(
treatment = as.factor(treatment),
genus = as.factor(genus),
wound_date = as.Date(as.character(wound_date), format = "%Y%m%d")
)
# Summary table
sample_info %>%
group_by(genus, treatment) %>%
summarize(n = n(), .groups = "drop") %>%
pivot_wider(names_from = treatment, values_from = n, names_prefix = "Treatment_") %>%
kable(caption = "Sample sizes by genus and treatment") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| genus | Treatment_0 | Treatment_1 | Treatment_2 |
|---|---|---|---|
| acr | 6 | 6 | 6 |
| por | 6 | 6 | 6 |
# Load processed weight data
initial_weight <- read_csv("data/processed/growth/initial_weight.csv", show_col_types = FALSE)
postwound_weight <- read_csv("data/processed/growth/postwound_weight.csv", show_col_types = FALSE)
day7_weight <- read_csv("data/processed/growth/day7postwound_weight.csv", show_col_types = FALSE)
final_weight <- read_csv("data/processed/growth/final_weight.csv", show_col_types = FALSE)
# Combine all weights
growth_data <- initial_weight %>%
select(coral_id, species, initial) %>%
left_join(postwound_weight %>% select(coral_id, species, postwound),
by = c("coral_id", "species")) %>%
left_join(day7_weight %>% select(coral_id, species, day7),
by = c("coral_id", "species")) %>%
left_join(final_weight %>% select(coral_id, species, final),
by = c("coral_id", "species"))
# Calculate growth rates
growth_data <- growth_data %>%
mutate(
growth_g = final - postwound,
growth_rate_g_day = growth_g / 23, # 23 days post-wounding
growth_normalized = growth_rate_g_day / initial
)
# Add treatment info
growth_data <- growth_data %>%
left_join(sample_info %>% select(coral_id, genus, treatment),
by = c("coral_id", "species" = "genus"))
# Summary
growth_data %>%
group_by(species, treatment) %>%
summarize(
mean_growth = mean(growth_rate_g_day, na.rm = TRUE),
se_growth = sd(growth_rate_g_day, na.rm = TRUE) / sqrt(n()),
.groups = "drop"
) %>%
kable(caption = "Growth rates by species and treatment (g/day)", digits = 4) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| species | treatment | mean_growth | se_growth |
|---|---|---|---|
| acr | 0 | 0.0743 | 0.0066 |
| acr | 1 | 0.0511 | 0.0087 |
| acr | 2 | 0.0516 | 0.0108 |
| por | 0 | 0.1910 | 0.0208 |
| por | 1 | 0.0872 | 0.0665 |
| por | 2 | 0.0535 | 0.0192 |
# Load averaged Fv/Fm data
pam_data <- read_csv("data/processed/pam/average_fvfm.csv", show_col_types = FALSE) %>%
select(-any_of("X")) %>% # Drop X column if it exists
mutate(
treatment = as.factor(treatment),
genus = as.factor(genus),
date = as.factor(date)
)
# Convert date to timepoint label
pam_data <- pam_data %>%
mutate(
timepoint = case_when(
date == "20230603" ~ "Day 7",
date == "20230619" ~ "Day 23",
TRUE ~ as.character(date)
),
timepoint = factor(timepoint, levels = c("Day 7", "Day 23"))
)
# Summary
pam_data %>%
group_by(genus, treatment, timepoint) %>%
summarize(
mean_fvfm = mean(average_fv_fm, na.rm = TRUE),
se_fvfm = sd(average_fv_fm, na.rm = TRUE) / sqrt(n()),
.groups = "drop"
) %>%
kable(caption = "Fv/Fm by genus, treatment, and timepoint", digits = 3) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| genus | treatment | timepoint | mean_fvfm | se_fvfm |
|---|---|---|---|---|
| acr | 0 | Day 7 | 697.222 | 10.119 |
| acr | 0 | Day 23 | 705.500 | 9.017 |
| acr | 1 | Day 7 | 702.000 | 8.572 |
| acr | 1 | Day 23 | 704.278 | 6.034 |
| acr | 2 | Day 7 | 699.333 | 9.146 |
| acr | 2 | Day 23 | 710.167 | 5.626 |
| por | 0 | Day 7 | 664.833 | 5.221 |
| por | 0 | Day 23 | 648.611 | 11.821 |
| por | 1 | Day 7 | 676.556 | 10.345 |
| por | 1 | Day 23 | 652.500 | 7.172 |
| por | 2 | Day 7 | 658.611 | 5.269 |
| por | 2 | Day 23 | 656.778 | 5.349 |
# Load surface area data
sa_data <- read_csv("data/processed/surface_area/final_surface_areas.csv", show_col_types = FALSE) %>%
select(-any_of(c("X", "...1"))) %>%
rename(coral_id = coral_number, genus_sa = taxa) %>%
mutate(genus_sa = tolower(genus_sa))
# Summary by genus
sa_data %>%
group_by(genus_sa) %>%
summarize(
mean_sa = mean(CSA_cm2, na.rm = TRUE),
sd_sa = sd(CSA_cm2, na.rm = TRUE),
min_sa = min(CSA_cm2, na.rm = TRUE),
max_sa = max(CSA_cm2, na.rm = TRUE),
.groups = "drop"
) %>%
kable(caption = "Surface area summary (cm²)", digits = 2) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| genus_sa | mean_sa | sd_sa | min_sa | max_sa |
|---|---|---|---|---|
| acropora | 17.88 | 5.10 | 11.55 | 30.27 |
| porites | 24.49 | 6.98 | 15.28 | 43.76 |
# Load respiration rate data for Porites
# Note: Acropora data needs to be processed first
# Function to load rates from a date
load_rates <- function(date, rate_type, genus = "Porites") {
path <- paste0("data/processed/respirometry/", genus, "/", rate_type, "/", date, "/")
if (!dir.exists(path)) {
return(NULL)
}
csv_files <- list.files(path, pattern = "\\.csv$", full.names = TRUE)
if (length(csv_files) == 0) {
return(NULL)
}
data <- read_csv(csv_files[1], show_col_types = FALSE) %>%
mutate(date = date, rate_type = rate_type) %>%
mutate(coral_id = as.numeric(coral_id)) # Standardize coral_id type
return(data)
}
# Load all Porites respiration data
resp_dates <- c("20230526", "20230528", "20230603", "20230619")
resp_data_list <- list()
for (date in resp_dates) {
resp <- load_rates(date, "Respiration", "Porites")
photo <- load_rates(date, "Photosynthesis", "Porites")
if (!is.null(resp)) resp_data_list[[paste0(date, "_resp")]] <- resp
if (!is.null(photo)) resp_data_list[[paste0(date, "_photo")]] <- photo
}
# Combine all data
if (length(resp_data_list) > 0) {
resp_all <- bind_rows(resp_data_list) %>%
clean_names() %>%
select(-starts_with("x")) %>%
filter(!is.na(coral_id))
# Convert date to timepoint
resp_all <- resp_all %>%
mutate(
timepoint = case_when(
date == "20230526" ~ "Pre-wound",
date == "20230528" ~ "Day 1",
date == "20230603" ~ "Day 7",
date == "20230619" ~ "Day 23",
TRUE ~ date
),
timepoint_num = case_when(
date == "20230526" ~ 0,
date == "20230528" ~ 1,
date == "20230603" ~ 7,
date == "20230619" ~ 23,
TRUE ~ as.numeric(NA)
),
timepoint = factor(timepoint, levels = c("Pre-wound", "Day 1", "Day 7", "Day 23"))
)
# Add genus and treatment info
resp_all <- resp_all %>%
mutate(genus = "por") %>%
left_join(sample_info %>% select(coral_id, treatment),
by = "coral_id")
cat("Loaded", nrow(resp_all), "respirometry measurements for Porites\n")
cat("Dates:", paste(unique(resp_all$date), collapse = ", "), "\n")
cat("Rate types:", paste(unique(resp_all$rate_type), collapse = ", "), "\n")
} else {
cat("No respirometry data found\n")
resp_all <- NULL
}## Loaded 304 respirometry measurements for Porites
## Dates: 20230526, 20230528, 20230603, 20230619
## Rate types: Respiration, Photosynthesis
# Growth by treatment and species
p1 <- ggplot(growth_data, aes(x = treatment, y = growth_rate_g_day, fill = treatment)) +
geom_boxplot(alpha = 0.7) +
geom_jitter(width = 0.2, alpha = 0.5, size = 2) +
facet_wrap(~species, labeller = labeller(species = c("acr" = "Acropora", "por" = "Porites"))) +
scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
labs(
title = "Growth Rate by Treatment and Species",
x = "Treatment",
y = "Growth Rate (g/day)"
) +
theme_few(base_size = 12)
# Normalized growth
p2 <- ggplot(growth_data, aes(x = treatment, y = growth_normalized, fill = treatment)) +
geom_boxplot(alpha = 0.7) +
geom_jitter(width = 0.2, alpha = 0.5, size = 2) +
facet_wrap(~species, labeller = labeller(species = c("acr" = "Acropora", "por" = "Porites"))) +
scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
labs(
title = "Size-Normalized Growth Rate",
x = "Treatment",
y = "Normalized Growth Rate (g/day / initial mass)"
) +
theme_few(base_size = 12)
print(p1)# Fv/Fm over time by treatment
p3 <- ggplot(pam_data, aes(x = timepoint, y = average_fv_fm, fill = treatment)) +
geom_boxplot(alpha = 0.7) +
geom_jitter(width = 0.2, alpha = 0.5, size = 2) +
facet_wrap(~genus, labeller = labeller(genus = c("acr" = "Acropora", "por" = "Porites"))) +
scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
labs(
title = "Photosynthetic Efficiency (Fv/Fm) Over Time",
x = "Timepoint",
y = "Fv/Fm"
) +
theme_few(base_size = 12)
# By treatment across species
p4 <- ggplot(pam_data, aes(x = genus, y = average_fv_fm, fill = treatment)) +
geom_boxplot(alpha = 0.7) +
facet_wrap(~timepoint) +
scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
labs(
title = "Fv/Fm by Species and Treatment",
x = "Genus",
y = "Fv/Fm"
) +
scale_x_discrete(labels = c("acr" = "Acropora", "por" = "Porites")) +
theme_few(base_size = 12)
print(p3)if (!is.null(resp_all)) {
# Filter out blank corals (coral_id == 0 or 1)
resp_filtered <- resp_all %>%
filter(!coral_id %in% c(0, 1))
# Respiration over time
resp_only <- resp_filtered %>% filter(rate_type == "Respiration")
p5 <- ggplot(resp_only, aes(x = timepoint, y = umol_l_sec, fill = treatment)) +
geom_boxplot(alpha = 0.7) +
geom_jitter(width = 0.2, alpha = 0.3, size = 1) +
scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
labs(
title = "Respiration Rate Over Time (Porites)",
subtitle = "Negative values indicate O2 consumption",
x = "Timepoint",
y = "Respiration Rate (µmol O2 / L / sec)"
) +
theme_few(base_size = 12)
# Photosynthesis over time
photo_only <- resp_filtered %>% filter(rate_type == "Photosynthesis")
p6 <- ggplot(photo_only, aes(x = timepoint, y = umol_l_sec, fill = treatment)) +
geom_boxplot(alpha = 0.7) +
geom_jitter(width = 0.2, alpha = 0.3, size = 1) +
scale_fill_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
labs(
title = "Photosynthesis Rate Over Time (Porites)",
subtitle = "Positive values indicate O2 production",
x = "Timepoint",
y = "Photosynthesis Rate (µmol O2 / L / sec)"
) +
theme_few(base_size = 12)
print(p5)
print(p6)
ggsave("reports/Figures/respiration_timeseries.png", p5, width = 10, height = 6, dpi = 300)
ggsave("reports/Figures/photosynthesis_timeseries.png", p6, width = 10, height = 6, dpi = 300)
# Combined P and R plot
p7 <- ggplot(resp_filtered, aes(x = timepoint_num, y = umol_l_sec, color = treatment)) +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.5, position = position_dodge(width = 1)) +
stat_summary(fun = "mean", geom = "point", size = 3, position = position_dodge(width = 1)) +
stat_summary(fun = "mean", geom = "line", position = position_dodge(width = 1), aes(group = treatment)) +
facet_wrap(~rate_type, scales = "free_y") +
scale_color_manual(values = treatment_colors, labels = treatment_labels, name = "Treatment") +
labs(
title = "O2 Exchange Rates Over Time (Porites)",
x = "Days Post-Wounding",
y = "Rate (µmol O2 / L / sec)"
) +
theme_few(base_size = 12)
print(p7)
ggsave("reports/Figures/O2_rates_combined.png", p7, width = 12, height = 6, dpi = 300)
}##
## === PORITES GROWTH ANALYSIS ===
por_growth <- growth_data %>% filter(species == "por")
# Linear model
mod_growth_por <- lm(growth_normalized ~ treatment, data = por_growth)
print(anova(mod_growth_por))## Analysis of Variance Table
##
## Response: growth_normalized
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 2 0.00010679 5.3393e-05 3.8337 0.0452 *
## Residuals 15 0.00020891 1.3927e-05
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call:
## lm(formula = growth_normalized ~ treatment, data = por_growth)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0027700 -0.0019869 -0.0007704 0.0002741 0.0120689
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.008133 0.001524 5.338 8.28e-05 ***
## treatment1 -0.004687 0.002155 -2.175 0.0460 *
## treatment2 -0.005540 0.002155 -2.571 0.0213 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003732 on 15 degrees of freedom
## Multiple R-squared: 0.3383, Adjusted R-squared: 0.25
## F-statistic: 3.834 on 2 and 15 DF, p-value: 0.0452
# Pairwise comparisons
emm_growth_por <- emmeans(mod_growth_por, ~treatment)
pairs_growth_por <- pairs(emm_growth_por, adjust = "tukey")
print(pairs_growth_por)## contrast estimate SE df t.ratio p.value
## treatment0 - treatment1 0.004687 0.00215 15 2.175 0.1082
## treatment0 - treatment2 0.005540 0.00215 15 2.571 0.0525
## treatment1 - treatment2 0.000853 0.00215 15 0.396 0.9177
##
## P value adjustment: tukey method for comparing a family of 3 estimates
##
## === ACROPORA GROWTH ANALYSIS ===
acr_growth <- growth_data %>% filter(species == "acr")
mod_growth_acr <- lm(growth_normalized ~ treatment, data = acr_growth)
print(anova(mod_growth_acr))## Analysis of Variance Table
##
## Response: growth_normalized
## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 2 6.0109e-05 3.0054e-05 4.0222 0.03994 *
## Residuals 15 1.1208e-04 7.4722e-06
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Call:
## lm(formula = growth_normalized ~ treatment, data = acr_growth)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0039210 -0.0025156 -0.0004061 0.0023428 0.0036627
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.010643 0.001116 9.537 9.29e-08 ***
## treatment1 -0.003605 0.001578 -2.285 0.0373 *
## treatment2 -0.004100 0.001578 -2.598 0.0202 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.002734 on 15 degrees of freedom
## Multiple R-squared: 0.3491, Adjusted R-squared: 0.2623
## F-statistic: 4.022 on 2 and 15 DF, p-value: 0.03994
emm_growth_acr <- emmeans(mod_growth_acr, ~treatment)
pairs_growth_acr <- pairs(emm_growth_acr, adjust = "tukey")
print(pairs_growth_acr)## contrast estimate SE df t.ratio p.value
## treatment0 - treatment1 0.003605 0.00158 15 2.285 0.0891
## treatment0 - treatment2 0.004100 0.00158 15 2.598 0.0500
## treatment1 - treatment2 0.000495 0.00158 15 0.313 0.9475
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# Linear mixed model with timepoint as fixed effect and coral as random effect
cat("\n=== PORITES FV/FM ANALYSIS ===\n")##
## === PORITES FV/FM ANALYSIS ===
por_pam <- pam_data %>% filter(genus == "por")
mod_pam_por <- lmer(average_fv_fm ~ treatment * timepoint + (1|coral_id), data = por_pam)
print(anova(mod_pam_por, type = 3))## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 206.01 103.00 2 15 0.4192 0.66503
## timepoint 1773.35 1773.35 1 15 7.2175 0.01691 *
## treatment:timepoint 762.23 381.11 2 15 1.5511 0.24417
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## timepoint = Day 7:
## treatment emmean SE df lower.CL upper.CL
## 0 665 7.98 26.6 648 681
## 1 677 7.98 26.6 660 693
## 2 659 7.98 26.6 642 675
##
## timepoint = Day 23:
## treatment emmean SE df lower.CL upper.CL
## 0 649 7.98 26.6 632 665
## 1 652 7.98 26.6 636 669
## 2 657 7.98 26.6 640 673
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## === ACROPORA FV/FM ANALYSIS ===
acr_pam <- pam_data %>% filter(genus == "acr")
mod_pam_acr <- lmer(average_fv_fm ~ treatment * timepoint + (1|coral_id), data = acr_pam)
print(anova(mod_pam_acr, type = 3))## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 34.45 17.23 2 15 0.0632 0.9390
## timepoint 457.48 457.48 1 15 1.6792 0.2146
## treatment:timepoint 115.73 57.86 2 15 0.2124 0.8110
## timepoint = Day 7:
## treatment emmean SE df lower.CL upper.CL
## 0 697 8.26 27 680 714
## 1 702 8.26 27 685 719
## 2 699 8.26 27 682 716
##
## timepoint = Day 23:
## treatment emmean SE df lower.CL upper.CL
## 0 706 8.26 27 689 722
## 1 704 8.26 27 687 721
## 2 710 8.26 27 693 727
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
if (!is.null(resp_all)) {
# Filter Porites, remove blanks, remove pre-wound timepoint
resp_analysis <- resp_filtered %>%
filter(timepoint != "Pre-wound")
# Separate analysis for Respiration
cat("\n=== PORITES RESPIRATION ANALYSIS ===\n")
resp_resp <- resp_analysis %>% filter(rate_type == "Respiration")
# Linear mixed model
mod_resp <- lmer(umol_l_sec ~ treatment * timepoint_num + (1|coral_id),
data = resp_resp,
contrasts = list(treatment = "contr.sum"))
print(anova(mod_resp, type = 3))
# Pairwise comparisons by timepoint
cat("\nPairwise comparisons at each timepoint:\n")
for (tp in unique(resp_resp$timepoint)) {
cat("\n", tp, ":\n")
resp_tp <- resp_resp %>% filter(timepoint == tp)
if (nrow(resp_tp) > 0) {
mod_tp <- lm(umol_l_sec ~ treatment, data = resp_tp)
emm_tp <- emmeans(mod_tp, ~treatment)
print(pairs(emm_tp, adjust = "tukey"))
}
}
# Analysis for Photosynthesis
cat("\n=== PORITES PHOTOSYNTHESIS ANALYSIS ===\n")
resp_photo <- resp_analysis %>% filter(rate_type == "Photosynthesis")
mod_photo <- lmer(umol_l_sec ~ treatment * timepoint_num + (1|coral_id),
data = resp_photo,
contrasts = list(treatment = "contr.sum"))
print(anova(mod_photo, type = 3))
# Pairwise comparisons
cat("\nPairwise comparisons at each timepoint:\n")
for (tp in unique(resp_photo$timepoint)) {
cat("\n", tp, ":\n")
photo_tp <- resp_photo %>% filter(timepoint == tp)
if (nrow(photo_tp) > 0) {
mod_tp <- lm(umol_l_sec ~ treatment, data = photo_tp)
emm_tp <- emmeans(mod_tp, ~treatment)
print(pairs(emm_tp, adjust = "tukey"))
}
}
}##
## === PORITES RESPIRATION ANALYSIS ===
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 1.5530e-05 7.7651e-06 2 97.386 0.7807 0.4609
## timepoint_num 1.7065e-05 1.7065e-05 1 85.905 1.7157 0.1937
## treatment:timepoint_num 9.9933e-06 4.9967e-06 2 85.905 0.5024 0.6069
##
## Pairwise comparisons at each timepoint:
##
## Day 1 :
## contrast estimate SE df t.ratio p.value
## treatment0 - treatment1 -0.000253 0.00037 33 -0.684 0.7742
## treatment0 - treatment2 0.000245 0.00037 33 0.661 0.7877
## treatment1 - treatment2 0.000498 0.00037 33 1.345 0.3811
##
## P value adjustment: tukey method for comparing a family of 3 estimates
##
## Day 7 :
## contrast estimate SE df t.ratio p.value
## treatment0 - treatment1 0.00277 0.00223 33 1.245 0.4357
## treatment0 - treatment2 -0.00121 0.00223 33 -0.542 0.8515
## treatment1 - treatment2 -0.00398 0.00223 33 -1.787 0.1896
##
## P value adjustment: tukey method for comparing a family of 3 estimates
##
## Day 23 :
## contrast estimate SE df t.ratio p.value
## treatment0 - treatment1 -0.000617 0.000322 33 -1.918 0.1495
## treatment0 - treatment2 0.000285 0.000322 33 0.886 0.6527
## treatment1 - treatment2 0.000903 0.000322 33 2.804 0.0222
##
## P value adjustment: tukey method for comparing a family of 3 estimates
##
## === PORITES PHOTOSYNTHESIS ANALYSIS ===
## Type III Analysis of Variance Table with Satterthwaite's method
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## treatment 4.0666e-05 2.0333e-05 2 96.201 1.5586 0.2157
## timepoint_num 2.7186e-05 2.7186e-05 1 83.987 2.0839 0.1526
## treatment:timepoint_num 8.7910e-06 4.3957e-06 2 83.987 0.3369 0.7149
##
## Pairwise comparisons at each timepoint:
##
## Day 1 :
## contrast estimate SE df t.ratio p.value
## treatment0 - treatment1 -4.15e-05 0.000245 33 -0.169 0.9843
## treatment0 - treatment2 2.83e-04 0.000245 33 1.156 0.4875
## treatment1 - treatment2 3.25e-04 0.000245 33 1.325 0.3917
##
## P value adjustment: tukey method for comparing a family of 3 estimates
##
## Day 7 :
## contrast estimate SE df t.ratio p.value
## treatment0 - treatment1 -0.004650 0.0024 33 -1.934 0.1450
## treatment0 - treatment2 0.000857 0.0024 33 0.357 0.9325
## treatment1 - treatment2 0.005506 0.0024 33 2.291 0.0712
##
## P value adjustment: tukey method for comparing a family of 3 estimates
##
## Day 23 :
## contrast estimate SE df t.ratio p.value
## treatment0 - treatment1 -0.000330 0.000367 33 -0.901 0.6436
## treatment0 - treatment2 -0.000168 0.000367 33 -0.458 0.8914
## treatment1 - treatment2 0.000163 0.000367 33 0.443 0.8977
##
## P value adjustment: tukey method for comparing a family of 3 estimates
# Overall sample summary
sample_info %>%
group_by(genus, treatment) %>%
summarize(n = n(), .groups = "drop") %>%
pivot_wider(names_from = genus, values_from = n) %>%
kable(caption = "Sample sizes by treatment and genus",
col.names = c("Treatment", "Acropora", "Porites")) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))| Treatment | Acropora | Porites |
|---|---|---|
| 0 | 6 | 6 |
| 1 | 6 | 6 |
| 2 | 6 | 6 |
if (!is.null(resp_all)) {
# Calculate mean rates by treatment, timepoint, and rate type
summary_stats <- resp_filtered %>%
group_by(treatment, timepoint, rate_type) %>%
summarize(
n = n(),
mean_rate = mean(umol_l_sec, na.rm = TRUE),
se_rate = sd(umol_l_sec, na.rm = TRUE) / sqrt(n()),
.groups = "drop"
)
summary_stats %>%
pivot_wider(
names_from = rate_type,
values_from = c(mean_rate, se_rate),
names_sep = "_"
) %>%
kable(caption = "Mean O2 exchange rates (µmol/L/sec) ± SE", digits = 4) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
scroll_box(width = "100%", height = "400px")
}| treatment | timepoint | n | mean_rate_Photosynthesis | mean_rate_Respiration | se_rate_Photosynthesis | se_rate_Respiration |
|---|---|---|---|---|---|---|
| 0 | Pre-wound | 12 | -0.0011 | -0.0031 | 0.0003 | 0.0005 |
| 0 | Day 1 | 12 | -0.0012 | -0.0033 | 0.0002 | 0.0003 |
| 0 | Day 7 | 12 | 0.0007 | -0.0037 | 0.0005 | 0.0013 |
| 0 | Day 23 | 12 | 0.0005 | -0.0028 | 0.0002 | 0.0003 |
| 1 | Pre-wound | 12 | -0.0001 | -0.0029 | 0.0003 | 0.0002 |
| 1 | Day 1 | 12 | -0.0012 | -0.0031 | 0.0002 | 0.0002 |
| 1 | Day 7 | 12 | 0.0053 | -0.0065 | 0.0029 | 0.0024 |
| 1 | Day 23 | 12 | 0.0008 | -0.0022 | 0.0002 | 0.0003 |
| 2 | Pre-wound | 12 | 0.0000 | -0.0029 | 0.0003 | 0.0002 |
| 2 | Day 1 | 12 | -0.0015 | -0.0036 | 0.0002 | 0.0003 |
| 2 | Day 7 | 12 | -0.0002 | -0.0025 | 0.0003 | 0.0002 |
| 2 | Day 23 | 12 | 0.0006 | -0.0031 | 0.0003 | 0.0002 |
Analysis completed: 2025-10-28 R version: R version 4.5.1 (2025-06-13)